Utilizing experimental design and desirability function in optimizing RP-HPLC method for simultaneous determination of some skeletal muscle relaxants and analgesics

An experimental design and response surface methodologies using Plackett–Burman and Box-Behnken designs were applied for selecting and optimizing the most appropriate parameters which significantly affect the separation and quantitative estimation of five skeletal muscle relaxants and four analgesic drugs (baclofen, methocarbamol, dantrolene sodium, orphenadrine citrate, cyclobenzaprine hydrochloride, ketoprofen, etoricoxib, ibuprofen, and mefenamic acid) with a relatively short duration of analysis in a single run. For the separation of the nine drugs, an INERTSIL ODS-V3-5 µm C18 column (250 × 4.6 mm I.D.) was used with the optimum mobile phase conditions (45.15 mM ammonium acetate buffer pH 5.56 adjusted with acetic acid, acetonitrile, and methanol in a ratio of 30.5:29.5:40, v/v/v with a flow rate of 1.5 mL/min) and UV-detection at 220 nm. The optimized method was successfully subjected to the validation steps as described in ICH guidelines for linearity, precision, accuracy, robustness, and sensitivity. The optimized and validated method was effectively applied to determine the content of the studied drugs in their pharmaceutical preparations and to expand its applicability to the counterfeit estimation of etoricoxib in different brands of tablet dosage forms.

The prevalence of counterfeit and substandard drugs is becoming a common worldwide problem which has been rapidly growing in recent years and affecting public health according to the WHO reports 1 .The term "counterfeit drug" refers to drugs that have been deliberately and fraudulently manufactured with mislabeling regarding their source and/or identity to become as genuine drugs such as products with the wrong or the right components, inadequate or insufficient amounts of the component(s), fraudulent packaging, or products without active ingredients 2 .The major causes of the prevalence of substandard and counterfeit drugs are the absence of government authorities' surveillance and monitoring; lack of awareness about the dangers of fake and fraudulent drugs, especially by pharmacists and patients; the online sale of counterfeit drugs; and the economic and political situation of the country, e.g.low economic situation due to conflicts and wars 3 .The planned treatment goals and outcomes can't be achieved if patients take counterfeit or substandard drugs.This could significantly impact their health and overall quality of life, potentially resulting in harm or even loss of life.Patients are not the only ones impacted by counterfeit and substandard pharmaceuticals; the entire health care system, medical professionals, the drug companies, drug authorities, Ministries of Health, and national economies are also affected.Common target drugs for counterfeiting include high-demand, pricey medications, e.g., chemotherapeutic agents, antivirals, antibiotics, erectile dysfunction medications, analgesics, antihistamines, weight loss aids, hormones, steroids, and antianxiety medications 4 .
Skeletal muscle relaxants have been proven to be highly effective therapeutic drugs for treating spasticity and muscle spasm, as well as being employed as surgical anesthetic adjuncts.The widespread use of skeletal muscle relaxants in medicine encouraged the creation (development) of simple, precise, accurate, sensitive, and easy-touse methods for determining their concentration in pure prepared mixtures and pharmaceutical dosage forms 5 .
Baclofen (BCL), chemically named as β-(aminomethyl)-p-chlorohydrocinnamic acid, is a selective agonist of the GABA-B (γ -aminobutyric acid) receptor because of its structural resemblance to the neurotransmitter γ-aminobutyric acid, which is used as a centrally acting muscle relaxant for the symptomatic treatment of a variety of neurological and musculoskeletal disorders, including muscle spasticity and multiple sclerosis 6 .Methocarbamol (MTL) is a centrally acting muscle relaxant with the chemical name ( ±)-3-(o-methoxyphenoxy)-1,2-propanediol-1-carbamate 7 , which is commonly used for the short-term symptomatic treatment of muscle spasms and is also used in combination with non-steroidal anti-inflammatory drugs such as ibuprofen to treat musculoskeletal disorders such as painful muscle spasms 8 .Orphenadrine citrate (ORF), chemically known as ( ±)-N,N-Dimethyl-2-[(o-methyl-α phenylbenzyl)oxyethylamine citrate, is a centrally acting skeletal muscle relaxant with anti-cholinergic (anti-muscarinic) and anti-histaminic effects which is used to alleviate painful muscle spasm 9 .Dantrolene sodium (DAN), chemically named 1-[ [5-(p-Nitrophenyl)furfurylidene]amino]hydantoin sodium salt hydrate, is a directly acting skeletal muscle relaxant which is used intravenously for treating and preventing malignant hyperthermia and related conditions, such as neuroleptic malignant syndrome, and also orally for controlling muscle spasticity 10,11 .Cyclobenzaprine hydrochloride (CYC), chemically designated as 3-(5H-dibenzo[a,d]cyclohepten-5-ylidene)-N,N-dimethyl-1-propanamine hydrochloride, is an efficient skeletal muscle relaxant with a central action.It works primarily in the brain stem to reduce somatic motor activity, which affects the alpha and gamma motor systems.Additionally, it is utilized to treat painful muscular spasms that are frequently linked to musculoskeletal disorders 12 .Ketoprofen (KPN), chemically designated as (RS)-2-(3benzoylphenyl)propionic acid, is a derivative of propionic acid with anti-inflammatory, analgesic, and antipyretic activities and possesses a non-selective inhibitory effect toward cyclooxygenase enzymes.It is used to treat mild to moderate pain in conditions like dysmenorrhea, migraines, post-operative pain, dental pain, and musculoskeletal and joint disorders like osteoarthritis and rheumatoid arthritis 13,14 .Etoricoxib (ETX), chemically known as 5-chloro-6'-methyl-3-[4-(methylsulfonyl)phenyl]-2,3'-bipyridine, is a non-steroidal anti-inflammatory drug (NSAID) that selectively inhibits the COX-2 enzyme and represents the second generation of COX-2 inhibitors.It is currently used to treat a number of painful disorders, including osteoarthritis, acute gout, ankylosing spondylitis, rheumatoid arthritis, dysmenorrhea, low back pain, short-term treatment of moderate pain after dental surgery, and acute postoperative pain in adults.When compared to nonselective NSAIDs, etoricoxib is typically more expensive 15,16 .Ibuprofen (IPN), chemically named (RS)-2-(4-(2-methylpropyl)phenyl)propanoic acid 17 , is a derivative of propionic acid that possesses a non-selective inhibitory effect toward cyclooxygenase enzymes and exhibits anti-inflammatory, analgesic, and antipyretic activities.It is licensed for the relief of pain and inflammation in rheumatic disease (juvenile idiopathic arthritis), mild to moderate dysmenorrhea, postoperative analgesia, migraine, dental pain, soft tissue injuries, pyrexia with discomfort, post-immunization pyrexia in infants, and other musculoskeletal disorders 13 .Mefenamic acid (MFC), chemically known as N-2,3-xylylanthranilic acid, is an anthranilic acid derivative with potent anti-inflammatory, antipyretic, and analgesic activities that is used for treating symptoms associated with menstruation, such as dysmenorrhea and menorrhagia, as well as pain and inflammation associated with osteoarthritis, rheumatoid arthritis, postoperative pain, and dental pain 18,19 .All the chemical structures of the cited drugs are shown in Fig. 1.
The development of chromatographic methods is influenced by various variables, and optimizing them is crucial for developing methods that are reliable, precise, and accurate.These methods should yield desired outcomes aligned with predefined objectives and the intended application.In accordance with ICH Q14 guidelines, there will be an increased emphasis on understanding and controlling the behavior of method variables during analytical method development.Applying the analytical quality by design (AQbD) approach, based on quality risk management (QRM) and design of experiments (DoE) principles, facilitates the understanding and control of potential method variables (PMV).Integrating DoE into chromatographic method development enhances comprehension of the relationships between critical method variables (CMVs) and critical analytical attributes (CAAs).It enables systematic exploration of the analytical design space (ADS), facilitating method optimization, robustness testing, and the identification of optimal conditions for achieving high-quality analytical outcomes.Optimizing CMVs is achieved through Response Surface Modeling (RSM) based on DoE, with exploration of the ADS ensuring the attainment of desired CAAs in accordance with the specified analytical target profile (ATP) for the chromatographic method 20 .
Plackett-Burman designs (PBD) are highly effective for screening experiments, but they are less adaptable to smoothly transitioning into optimization phases compared to two-level factorial designs.Nonetheless, their efficiency in assessing main effects with a minimal number of experimental runs makes them well-suited for robustness studies, especially when exploring slight variations from method conditions 21 .
Box-Behnken designs (BBD), a specific type of response surface methodology (RSM), demonstrate efficiency by requiring fewer experimental runs compared to designs like the central composite design with the same number of factors.This cost-effectiveness enables the effective estimation of both first-and second-order coefficients in the model, providing comprehensive insights into the relationships between factors and the response(s).However, unlike certain other RSM designs, they do not directly incorporate runs from a factorial experiment and necessitate a minimum of three factors for proper implementation 22 .
The desirability function enables the analyst to assess how effectively a combination of experimental conditions (factor levels) meets the predefined objectives for the responses.This involves achieving the optimal value for all assessed responses simultaneously and considering the researcher's priorities during the optimization process.The individual desirability (d) assesses how the settings optimize a single response, while the composite desirability (D) assesses how the settings optimize a set of responses overall.Desirability ranges from zero to one, with 1 representing the ideal combination of response values where all are within acceptable limits and zero indicating that one or more responses fall outside acceptable limits 23 .
There have been no attempts to provide a single method for separating and determining a variety of skeletal muscle relaxants with some NSAIDs.Our goal is to develop a new, versatile, economic, simple, rapid, and robust RP-HPLC method for separation and qualitative and quantitative analysis of the nine drugs in a single process.This will be achieved by employing a DoE approach and RSM, specifically utilizing PBD for screening and BBD for optimizing and estimating the robustness of the method.Statistical regression analysis was applied to investigate how different chromatographic conditions affected the separation of these drugs.The developed method was also designed with the goal of being employed in the detection and determination of any massage creams counterfeited with skeletal muscle relaxants, as well as for the counterfeit detection of etoricoxib in different brands of tablet dosage forms sold in Yemen.It is advised to use it in the quality control check for potential counterfeits of these drugs.We advise pharmaceutical companies or drug authorities to use a single analytical method to verify the quality of various pharmaceutical products during post-market monitoring because there is a huge number of different products that are taken into account for evaluation in these situations 58 .

Experimental Instrumentation
The HPLC system used was an Agilent-1100 series from Japan with a G1310A-ISO pump and a G1314A-variable wavelength detector.For data collection, analysis, and instrument control, Agilent Chemstation was used.A C18 column (INERTSIL ODS-3 V, 250 mm × 4.6, 5 µm) was used for HPLC separation.A sonicator was used (Soniclean 120 T, Barton, SA, Australia).The pH was measured and adjusted using a pH meter (Jenway-3505, Essex, UK).For sample preparation and solvent filtration, non-sterile Minisart NY syringe filters and membrane filters (0.22 m and 0.45 m, SARTORIUS STEDIM BIOTECH, Germany) were used.PURELAB Option-Q was used to obtain ultra-pure water (ELGA LABWATER, UK).Data optimization and DoE were done using the Minitab-17 (Minitab, Inc., State College, Pennsylvania, USA).

Preparation of stock and working samples solutions
Tablets dosage forms sample solutions Ten tablets of MYLOBAC, METHOCARBEX, ARCOXIA, or MULTI-RELAX were individually weighed and finely pulverized to form a homogenous powder.An accurately weighted quantity of the powder equivalent to one tablet of each dosage form was separately transferred into a 100-mL volumetric flask.To ensure that the drugs were completely dissolved, 50 mL of mobile phase was added to each flask and sonicated for at least 15 min.After cooling, each flask was completed to the mark with the mobile phase.The final step was to filter the resulting solutions through Whatman filter papers, discard the first few milliliters, and then filter them again through 0.22 μm syringe filters.

Capsules dosage forms sample solutions
The contents of 10 capsules of each of KETOLGIN, DANTRELAX, or MEFENTAN were weighed and completely mixed.An accurate portion, equivalent to the content of one capsule, was placed in a 100-mL volumetric flask.To ensure that the drugs were completely dissolved, 50 mL of mobile phase was added to each flask and shaken well by using a sonicator for at least 15 min.After cooling, each solution was diluted and completed to volume with the mobile phase.Finally, the resulting solutions were passed through Whatman filter papers, and the first few mLs were discarded and then filtered again using 0.22-μm syringe filters.

Ampoules sample solutions
Ten NORFLEX ampoules were thoroughly mixed.An equivalent amount of one ampoule was accurately introduced into a 100-mL volumetric flask.The volume was adjusted to the mark with the mobile phase.

Massage creams sample solutions
Five massage cream samples were bought from the Yemeni and Egyptian markets.One gram from each of the five massage cream products was introduced separately into a 50-mL volumetric flask.Then, 30 mL of mobile phase was added, and the solution was mixed with a vortex mixer and sonicated for 15 min to disperse the sample thoroughly.The volume was then completed with the mobile phase to the mark and passed through a suitable filter of 0.45 µm filter paper and a 2-µm pore size using a syringe filter.The first few mLs were discarded.An aliquot of 2 ml from each solution was introduced into a set of 10 ml volumetric flasks, and it was then completed to volume with the mobile phase.

Tablets dosage forms sample solutions for counterfeit verification
Ten tablets of ETRICIB-90, ETROPAIN-90, ENTORIK-90, ETOHEAL-90, ARCORAR-90, or E-COX-90 were weighed individually and finely pulverized to form a homogenous powder.An accurately weighted quantity of powder equivalent to one tablet of each dosage form was transferred into a 100-mL volumetric flask.To ensure that the drugs were completely dissolved, 50 mL of mobile phase was added to each flask and sonicated for at least 15 min.After cooling, each flask was completed to the mark with the mobile phase.Finally, the resulting solutions were filtered through Whatman filter papers, and the first few milliliters were discarded and then filtered again using 0.22 μm syringe filters.

Plackett-Burman design for screening
The PBD's objective was to determine, with the fewest possible runs, the significant effects of variables on the responses and details of these variables and responses are shown in Table 1.Initially, the main effects of four parameters were estimated using the 12 runs of PBD (Table 2), which also provided a measurement of the process's inherent variability and stability.

Box-Behnken design for factor optimization
After applying the PBD screening to select the most significant predictors on the responses (retention time of the last eluted drug and peak resolutions), the three-level four-factor BBD (Table 3), which has 27 runs with three center points, was used to assess the main and quadratic effects of the significant factors (the percentage of acetonitrile, pH of buffer, flow rate of mobile phase, and concentration of buffer) and their interaction effects on the chosen responses for further optimization of the HPLC analysis conditions.

Analysis of the cited drugs in their pharmaceutical dosage forms
For determining the nine studied drugs in their pharmaceutical dosage forms, the stock sample solutions were diluted with the mobile phase and then injected into the column in triplicate.To calculate the amounts of the cited drugs, the regression equations for each drug were used.

Analysis of the massage creams for counterfeit detection
After completing the extraction process and preparing the stock and working solutions of massage creams, 20 μl of each solution was injected into the column under the same ideal conditions that were chosen for the method, and then the resulting chromatogram was compared with the chromatograms resulting from the mixture of the cited drugs.

Analysis of the different brands of etoricoxib for counterfeit detection
For determining the counterfeit of ETX in six different brands of tablets, the stock sample solutions were diluted with the mobile phase and then injected into the column in triplicate.To calculate the amount of ETX by the standard addition method, the regression equation for ETX was used.

Test of system suitability
For determination of the system suitability of the proposed method, a standard mixture solution containing 50 µg/mL of the nine studied drugs was injected six times consecutively.Several parameters were examined, including column efficiency (theoretical plates), capacity factor, resolution, tailing factor, selectivity, and percentage relative standard deviation of retention times.

Method development using analytical quality by design (AQbD) approach
Chromatographic method development is impacted by a variety of method variables and optimizing them is essential for developing methods that are robust, accurate, and precise.These methods should produce desired outcomes in line with predefined objectives and the intended application of the method.In accordance with the ICH Q14 guidelines, there will be an increased focus on understanding and controlling the behavior of method variables during the development phase of an analytical method.The understanding and control of the potential method variable (PMV) can be achieved by applying the analytical quality by design (AQbD) approach, which is based on the principles of quality risk management (QRM) and design of experiments (DoE).The application of DoE in chromatographic method development enhances the understanding of the relationships between the critical method variables (CMVs) and the critical analytical attributes (CAAs).It facilitates a systematic exploration of the analytical design space (ADS), enabling method optimization, robustness testing, and the identification of optimal conditions for achieving high-quality analytical results 20 .The optimization of CMVs can be accomplished by employing Response Surface Modeling (RSM) based on DoE.Utilizing this approach to explore the ADS enables the achievement of desired CAAs in accordance with the specified Analytical Target Profile (ATP) for the chromatographic method.
In the initial phase of implementing AQbD for the RP-HPLC method, the focus was on precisely defining the analytical target profile (ATP) and critical analytical attributes (CAAs).
The method's analytical target profile focuses on developing and validating a versatile and robust RP-HPLC method for simultaneously estimating five skeletal muscle relaxants and four analgesic drugs.The objective is to achieve a good resolution between chromatographic peaks exceeding 1.5 with a minimum run time.Consequently, the Critical Analytical Attributes (CAAs) selected for the development of the target RP-HPLC method include specific values for resolution between the studied drugs (R-1 to R-8) and retention time of the last eluted drug (t R-9) .

Experimental design for optimization of chromatographic conditions
Experimental design is sequentially used to plan the experiment runs and analyze the data in order to identify the significant variables (critical method variables) that need to be reconstructed and the factor settings that produce the best possible response(s) (critical analytical attributes).The main advantage of experimental design is that it requires less time, effort, and money because only a few experimental runs are carried out.The two steps of DoE have been applied: the screening step using BPD to determine which factors are important for explaining process variation; and the second step is the optimization step using response surface design, such as BBD to create a response surface that describes the relation between the responses (CAAs) and the variables (CMVs) and determines the optimum chromatographic conditions that provide a good resolution, symmetrical peak shape, and adequate retention time.

Screening variables using PBD
The primary objective of our method development is to efficiently identify the most influential factors among a set of potential method variables affecting the chromatographic separation process, i.e., to efficiently select critical factors that significantly impact the method's performance.Therefore, the PBD is well-suited for this initial screening phase because it allows us to assess the main effects of multiple factors with a relatively small number of experimental runs.This screening approach enables us to focus subsequent optimization efforts on the most influential factors, simplifying the overall method development process.In summary, the rationale for choosing the PBD lies in its effectiveness as a screening tool, allowing us to efficiently identify key factors for further investigation and optimization in the development of our chromatographic method.
In this study, the PBD was applied to examine the major impacts of four parameters (A, B, C, and D), on nine responses, which are presented in Table 1.These responses included t R-9 (retention time of MFC) and various resolution factors between peaks of the studied drugs.
The screening variables were set at two levels (− 1, 1) as depicted in Table 2. Twelve experimental runs, resulting from the combination of four variables, were designed and randomly conducted.Response data were collected based on the designated run order, as detailed in Table 2.
The significant main effects of the variables were statistically determined through analysis of variance (ANOVA) at a 95% confidence interval, as presented in the Supplementary Table S2.A graphical representation of the main effects was done using the main effect plots shown in Fig. 2 and the Pareto chart of standardized effects.In the Pareto chart, variables that are statistically significant are those that cross the reference line, as illustrated in Fig. 3.
Based on the ANOVA results and the Pareto chart (as shown in Fig. 3), A (ACN%), B (flow rate), and C (buffer pH) all had a significantly negative effect on tR-9 (retention time).R-1 was notably affected only by A (ACN%) and B (flow rate), while R-2 and R-3 were strongly influenced by all factors.R-4 was notably affected only by B (flow rate) and C (buffer pH), while R-5 and R-7 were significantly influenced by A (ACN%) and C (buffer pH).R-6 was significantly affected by B (flow rate), C (buffer pH), and D (acetate buffer concentration), whereas R-8 was significantly affected by A (ACN%), B (flow rate), and D (acetate buffer concentration).
According to the main effect plots (Fig. 2), the Pareto chart (Fig. 3), and the ANOVA results of the PBD (Supplementary file, Table S2), it was determined that all the studied factors (A, B, C, D) were significant.Therefore, they should be optimized in the BBD to explore their main and quadratic effects, along with their interaction effects, on the desired responses 59 .www.nature.com/scientificreports/Optimization of factors using the BBD BBD is a particular type of response surface methodology (RSM) that offers efficiency through a reduced number of experimental runs compared to other designs like the central composite design with the same number of factors.This makes them more cost-effective.They facilitate the effective estimation of both first-and secondorder coefficients in the model, providing deeper insights into the relationships between factors and the response.However, unlike some other RSM designs, they don't directly incorporate runs from a factorial experiment and require a minimum of three factors for proper implementation 22 .
In the response surface methodology, to extract a mathematical model depicting the relationship between factors and response variables, first-order and second-order response surface models are frequently employed.These models typically take the form of linear or quadratic polynomial functions for approximation.In simpler terms, the primary form of the first-order response surface model looks like this: Alternatively, the fundamental form of the second-order response surface model is expressed as follows: In this context, y represents the dependent variable (response); xi represents the ith component of the n-dimensional independent variable x; β 0 , βi, and βij denote the unknown parameters forming the column vector β; and ε represents the error term 60 .
The three-level, four-factor BBD with 27 runs, including three center points, was employed to optimize the experimental data and identify the ideal settings of HPLC parameters that result in good resolutions between the tested drugs with the shortest run time.The four factors selected for the optimization step were A (ACN%), B (flow rate), C (buffer pH), and D (acetate buffer concentration), and their levels are shown in Table 3.All experimental runs were carried out in a randomized order with three replications at the center points to assess the experimental error, reducing the effects of independent predictors (extraneous or uncontrollable conditions) that are not included in the study and may bias the observed results 59,61 .The inherent variance in parameters can also be estimated through randomization in order to properly draw statistical conclusions from the results of the experiment.The nine measured responses were t R-9 , and the eight resolutions (R-1, R-2, R-3, R-4, R-5, R-6, R-7, R-8) were between the nine drugs.As a result, nine regression models were built and used to determine the surface relationship between one or more factors and the response variables.The experimental runs and measured responses are shown in Table 3, which also provides a detailed description of the BBD matrix for method optimization of four independent factors and nine responses.After collecting data from each experimental point in a chosen design, it's essential to establish a mathematical equation that describes how the response behaves with respect to the different levels of the variables studied.This means we need to estimate the parameters (β) in Eqs. ( 1) and (2).The regression coefficients are calculated using the least squares method, which involves minimizing the sum of squared differences between the observed values and the predicted values.Subsequently, based on these regression coefficients, a response surface models are constructed, aiding in predicting the relationship between the variables in the dataset 62 .Using regression analysis for BBD, the model regression equations that describe the model as either linear or quadratic were derived and represented in Eqs.(3-11) as follows: where A, B, C and D are ACN%, flow rate, buffer pH, and acetate buffer concentration, respectively.AB, AC, AD, and CD indicate the interaction between the factors, while AA, BB, CC, and DD represent the quadratic terms for each factor.

Responses
The significant linear and quadratic effects of the studied factors on various measured responses, as well as their interactive effects, are summarized in Table 4.The positive sign of a coefficient in a regression equation indicates a positive correlation between the factor and the response, while a negative sign represents a negative correlation between the factor and the response.Regression equations and ANOVA analysis revealed significant models for multiple responses.For example, according to the data in Table 4, it is evident that tR-9, R-2, and R-3 have p-values lower than 0.05 for both linear, quadratic, and interaction models (the p-value less than 0.05 indicates that the effect is significant).Additionally, R-1 shows p-values lower than 0.05 for both linear and quadratic models, indicating significance for both.Also, the linear effect of the C factor is non-significant, but its interaction effect with A (AC) is significant.Meanwhile, R-8 exhibits p-values lower than 0.05 for both linear and interaction models.However, the selected models (equations) for all responses include main effects, interaction terms, and quadratic terms.This complexity arises because when multiple significant models exist for the same response (all with p-values less than 0.05), the selected model incorporates all significant terms from each model, resulting in a more complex model.The ANOVA test was used to statistically analyze the resulting regression models.The results show that the models for the selected independent variables are significant with P-values < 0.0001, and the results of the ANOVA analysis are shown in Table 4.The p-values for all models' lack of fit suggest that the data is well-fitted by these models (a p-value of more than 0.05 is required for the models' lack of fit), and the small p-values for the linear terms, the quadratic terms, and the interactions indicate that these effects are statistically significant (for the model and the variables to be significant, the p-value must be lower than 0.05).The small p-values for the interactions (p-value < 0.05) and the squared terms (p-value < 0.05) suggest there is curvature in the response surface.The high F-values for the models and factors suggest that the models are significant at a 95% level of confidence and that the models fit the data reasonably well.Backward elimination was performed using Minitab.In this process, Minitab initially includes all variables in the model and iteratively eliminates the non-significant variable at each step.The procedure continues until all variables in the model exhibit p-values less than or equal to the specified alpha-to-remove value (0.05).Table 4 displays the R-squared, adjusted R-squared, and predicted R-squared for each model.All models exhibited R-squared values ranging from 90 to 100%, indicating that they adequately explain 90% to 100% of the variances in responses and suggesting a relatively good fit to the data.The adjusted R-squared values ranged from 89 to 100%, and the predicted R-squared values ranged from 87 to 100%.The Adjusted R-squared value serves as a more reliable measure of the model's goodness of fit than the R-squared value, especially in situations where the model incorporates an excessive number of independent variables.High values of R-squared and adjusted R-squared indicate that the data have been fitted properly by the models 39 .As noted in Table 4, the predicted R-squared values were within accepted limits with regards to the adjusted R-squared values, i.e. the difference was less than 0.2, indicating that the experimental data and the fitted model are in good agreement and demonstrating the high predictive power of models (predicted R-squared should be as near to adjusted R-squared as possible in order to make precise predictions) 59,63 .
Examining residual plots was another method used to assess the model's fit adequacy and determine whether there are any outliers, non-normality, or non-random variation 61 .The residuals are normally distributed because all points form a straight line in the normal probability plot of the residuals, as illustrated in Fig. 4, and they are randomly scattered on either side of zero, demonstrating that the errors had a constant variance and were uncorrelated with one another, as shown in residual plots versus fits and versus order in Fig. 4.
The response surface was visually represented using 2D-contour plots (as shown in Fig. 5A and B).These plots illustrate how the response is influenced by each of the two variables and describe the shape of the response surface as shaded areas, contour lines, or both 64 and illustrate the impact of the interaction between pairs of factors on the response while keeping the other variables at their central point values.The curvature observed in these plots indicates the non-linear effects of the factors on the response.These graphical representations facilitate the prediction of the response at any point within the analytical design space (experimental field).The combined effect of A (ACN%) and B (flow rate) on responses (tR-9, R-1, R-2, R-3, R-4, R-5, R-7, and R-8) is illustrated in Fig. 5A,a1-a8.It is observed that an increase in ACN% from 20 to 30% leads to a decrease in the previous responses, except for R4, which increases.Additionally, an increase in the flow rate from 1 to 1.5 leads to a decrease in responses.The combined effect of A (ACN%) and C (buffer pH) on responses (t R-9 , R-2, R-3, R-4, R-5, R-7, and R-8) is illustrated in Fig. 5A,b1-b7.It is observed that an increase in buffer pH from 5.3 to 6.0 leads to a decrease in the responses t R-9 , R-2, R-5, and R-7, while it leads to an increase in R-3, R-4, and R-8.The combined effect of A (ACN%) and D (acetate buffer concentration) on responses (t R-9 , R-1, R-2, R-3, R-5, R-7, and R-8) is illustrated in Fig. 5A,c1-c8.In the case of the D (acetate buffer concentration), increasing the concentration from 10 to 50 mM reduces the previous responses (t R-9 , R2, R-5, R-7, and R-8) except for R-1 and R-3.In the beginning, it leads to an increase in these two responses to a certain point, and after that, the opposite happens, i.e., an increase in the buffer concentration, it leads to a decrease in the responses R-1 and R-3.
In order to calculate the composite or combined desirability (D), the individual desirabilities for each response must first be calculated.In general, a quick analysis time with good peak resolution is recommended so that the last optimization step is based on minimizing t R-9 , R-1, R-2, and R-8 and determining the remaining responses at target values as shown in Table 5.The following ideal conditions (the optimal solution) were optimized based on the satisfactory composite desirability that was obtained (D = 0.9998): 29.5 ACN%, flow rate of 1.5 mL/min, pH 5.56, and 45.15 mM acetate buffer, as can be seen in the optimization plot in the supplementary file, Fig. S1.For additional model validation and verification of the model's predictive ability, an experimental run using the ideal settings of the variables (the ideal solution) was carried out, with the observed results being compared to the predicted results, as shown in Supplementary Table S3.

Application of the optimized method
Application to the estimation the drugs' contents in laboratory prepared mixtures and pharmaceutical dosage forms The method was run under ideal chromatographic conditions for quantitative determination of the nine studied drugs in their pure standard mixtures and their pharmaceutical dosage forms, and all results were in good agreement with the expected values, as shown in Table 6.

Application to detection of counterfeits in massage creams
The optimized HPLC method was applied to detect whether the massage creams were adulterated with skeletal muscle relaxants and analgesic drugs by comparing the retention times of the peaks that were eluted from the massage creams to the retention times of the studied drug standards, as shown in Fig. 6.It was found that all massage cream products had no peaks at the retention times of the studied drug standards except ROTMOOV massage cream, which had a small peak at the retention time of CYC.When the studied ROTMOOV massage   www.nature.com/scientificreports/

Application to detection of counterfeits in etoricoxib tablets
The optimized method was applied to estimate the drugs' contents in 6 different brands of etoricoxib.According to the quantitative analysis of six different brands of etoricoxib coded (ETX-I to ETX-6) in Table 7, only ETX-6 (86.33% ± 0.76) had less than 90% of the active ingredient.In this situation, a drug's efficacy could be affected by an insufficient amount of the active ingredient.According to Indian pharmacopoeia, the percentage content of etoricoxib should fall within the range of 90%-110%.

System suitability tests
To determine the system suitability of the proposed method, several parameters were examined, with the results shown in Table 8.These parameters included column efficiency (theoretical plates), capacity factor, resolution, tailing factor, selectivity, and percentage relative standard deviation of retention times.www.nature.com/scientificreports/

Validation of the optimized method
The principles of the ICH recommendations have been applied to the validation of the proposed method 65 .

Linearity
Eight different concentrations of the nine studied drugs were analyzed in triplicates to achieve linearity.The good linearities were confirmed by high values for the regression coefficients.Nine calibration curves with nine regression equations were constructed by plotting the relationship between drug concentrations and the corresponding peak areas.All data from calibration curves, such as values of slopes, intercepts, and regression coefficients, as well as standard deviations for the intercept and slope, are shown in Table 6.

Accuracy
To validate the accuracy of the proposed method, five different concentrations of the nine studied drugs were injected three times.Following this, the percent recoveries of the studied drugs were calculated in their pure standard mixtures.The content recoveries of the studied drugs in formulated dosage forms were calculated by the standard addition technique.As seen in Table 6, good recovery values were obtained for drugs in both pure standard mixtures and in formulated dosage forms.

Precision
The intraday precision of the proposed method was estimated by analysis of three concentrations of

Detection and quantitation limits
Using the method based on the response standard deviation (Sa) and slope of the calibration curve (b), the detection and quantitation limits were calculated by using the following equations: LOD = 3.3 Sa/b and LOQ = 10 Sa/b.The good results reflect the good sensitivity of the proposed method as shown in Table 6.

Selectivity
The term "selectivity" refers to a method's capacity to detect and measure an analyte's response without interference from other analytes and excipients.It was evaluated by the examination of several laboratory-made mixtures of the nine studied drugs, as well as through the analysis of various dosage forms of the studied drugs.It is found that the drug peaks in the sample solutions of all dosage forms are similar to those in the standard solutions.Additionally, no chromatographic impurity interference was seen in any of the dosage forms.The good resolution between the drug peaks and the successful determination of the nine studied drugs in pharmaceutical dosage forms, as shown in Fig. 6 and Table 6, both confirmed the good selectivity.

Statistical analysis
The results obtained from the proposed method have been statistically analyzed using the Student's t-test and variance ratio F-test, at P = 0.05, with those obtained using the U.S.P. reference methods for BCL, MTL, KPN, DAN, ORF, CYC, IPN, and MFC 7 and the I.P reference method for ETX 66 .The results show no significant difference in accuracy and precision between the methods for each drug, as shown in Table 9.

Conclusion
In this work, screening, optimization, and robustness assessment were successfully carried out using experimental design (Plackett-Burman and Box-Behnken designs) to develop the RP HPLC method for the separation and simultaneous quantification of five drugs of skeletal muscle relaxants with four analgesic drugs (BCL, MTL, KPN, DAN, ORF, ETX, CYC, IPN, and MFC) in their pure standard mixtures and in pharmaceutical products.
The suggested method was quick, precise, accurate, and sensitive, and it allowed the separation of the nine drugs in less than 8 min of satisfactory runtime.The suggested method was successfully used to estimate the concentrations of the analyzed drugs in pharmaceutical products after being fully validated by the validation process described in the ICH recommendations.Moreover, the proposed method revealed a lower quantity of the active ingredient ETX in Indian pharmaceutical products than declared on the label.It also confirmed the absence of counterfeit drugs in the analyzed massage cream.The suggested method can therefore be applied, particularly in laboratories of quality control, for routine evaluation of pharmaceutical products containing the cited drugs and can also be applied to detect the potential counterfeiting of these drugs.

Figure 1 .
Figure 1.Chemical structures of the studied drugs.

Figure 2 .
Figure 2. Screening of the main effects of the studied factors on the measured responses.

Figure 3 .
Figure 3. Pareto chart of the standardized effects of the studied factors on the measured responses.

Figure 4 .
Figure 4. Normal probability plot, histogram, residuals versus fits and residuals versus order plots for the measured responses.

Figure 5 .Figure 5 .
Figure 5. (A) and (B) Contour plots showing the interaction effects of studied factors on the measured responses.

Table 1 .
The abbreviation and definition of the studied factors and responses.

Table 4 .
Reduced response models and statistical parameters obtained from ANOVA (after backward elimination) (AB, AC, AD, and CD indicate the interaction between the factors; AA, BB, CC, and DD represent the quadratic terms for each factor; F: F-value; P: P-value; The shaded rows represent final selected model).

Table 5 .
Criteria for the optimization of the responses.

Table 6 .
Results of the regression equations and validation parameters of the optimized analysis method (R 2 : regression coefficient; Sa standard deviation of intercept, Sb standard deviation of Slope, LOD limit of detection, LOQ limit of quantitation).

Table 7 .
Results of analysis of 6 brands of etoricoxib by the proposed method.

Table 8 .
Results of system suitability tests for the optimized HPLC method.